Critical behavior of the three-dimensional bond-diluted Ising spin glass: 
Finite-size scaling functions and universality 
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We study the three-dimensional (3D) bond-diluted Edwards- Anderson (EA) model with binary 
interactions at a bond occupation of 45% by Monte Carlo (MC) simulations. Using an efficient 
cluster MC algorithm we are able to determine the universal finite-size scaling (FSS) functions and 
the critical exponents with high statistical accuracy. We observe small corrections to scaling for 
the measured observables. The critical quantities and the FSS functions indicate clearly that the 
bond-diluted model for dilutions above the critical dilution p* , at which a spin glass (SG) phase 
appears, lies in the same universality class as the 3D undiluted EA model with binary interactions. 
A comparison with the FSS functions of the 3D site-diluted EA model with Gaussian interactions 
at a site occupation of 62.5% gives very strong evidence for the universality of the SG transition in 
the 3D EA model. 

PACS numbers: 75.10.Nr 
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I. INTRODUCTION 

Almost all the Monte Carlo (MC) simulations of the 
three-dimensional (3D) Edwards- Anderson (EA) model 1 
have been performed using either an undiluted Gaus- 
sian or binary coupling distribution.-- There are, how- 
ever, many reasons to study the 3D bond-diluted EA 
model with binary interactions. The simplest reason be- 
ing that it allows one, although in a very crude way, to 
approximate the presence of weak bonds in the model 
with Gaussian couplings without loosing the technical 
benefits of the binary couplings. The introduction of va- 
cant bonds reduces the frustration in the model, making 
it easier to find ground states, which should facilitate 
low-temperature MC simulations. Moreover, the bond- 
diluted model in a certain range of the dilution allows for 
an efficient use of replica cluster MC algorithms.— &£iiLI. 
It has been shown that the high-temperature series ex- 
pansion typically converges better for some intermediate 
dilution than for the undiluted case. 8 Similarly, it was 
found that the scaling of the ground state defect energies 
is much improved at intermediate dilutions in compari- 
son to the standard undiluted case, that showed to scale 
rather poorly. 9 This gives some evidence that in general 
the scaling violations of the diluted model, for some in- 
termediate range of the dilution, might be smaller than 
in the undiluted case. This fact is very important for 
spin glass (SG) simulations since the range of system 
sizes that are accessible to MC calculations is very limited 
due to the presence of large free energy barriers in the 
SG phase, causing slow thermalization. The problems of 
the very restricted simulation range are worsened by the 
fact that the leading correction-to-scaling exponent lo has 
been measured to be probably smaller than one , 10 ! 1 such 
that the leading corrections fall off only slowly. All this 
combined with the need of excellent statistics makes a re- 
liable estimate of critical exponents very difficult. In ad- 
dition to these rather technical aspects the bond-diluted 



± J EA model may be thought of simply as an EA model 
with a particular distribution of the couplings and there- 
fore it also serves to provide insights on the debated issue 
of the universality of the 3D EA model&12J^±2J&lL12. 

In this paper we first consider in detail the bond- 
diluted 3D EA model at a bond occupation of 45% 
which is close to the value where optimal scaling for the 
ground state defect energies was found* 9 - We first discuss 
the replica cluster MC algorithm we use and its perfor- 
mance, then we focus on the SG transition, the corre- 
sponding critical exponents and the universal finite-size 
scaling (FSS) functions of the correlation length £, the 
SG susceptibility xsg an d the Binder ratio g. We de- 
termine these FSS functions with very high statistical 
accuracy (up to 10 5 samples per lattice size) and use 
them to determine critical exponents with two different 
methods, namely, the quotient metho d 10 ! 19 and an ex- 
trapolation to infinite volume method^ We find excel- 
lent agreement between the critical exponents extracted 
by these two different methods, as well as with the val- 
ues given by the high-precision study of the undiluted 
binary model of Ballesteros et al*£ In addition, we find 
that the FSS functions agree very well with the ones ob- 
tained by Palassini and Caracciolo for the undiluted bi- 
nary modeliii This together gives strong evidence that 
the SG transition of the 3D bond-diluted EA model with 
binary couplings falls into the same universality class as 
the one of the corresponding undiluted model, as indi- 
cated by high-temperature expansion^ and as one might 
expect on general grounds^ We then compare the FSS 
functions to the ones of the 3D site-diluted EA model 
at a site occupation of 62.5% with a Gaussian distribu- 
tion of the couplings^ In contrast to the bond-diluted 
±J model the site-diluted Gaussian model shows very 
noticeable finite-size effects and therefore the agreement 
between the FSS functions of the two models emerges 
only at large enough sizes. Despite the fact that the two 
models have finite-size effects of clearly different magni- 
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tude the agreement between the FSS functions gets ex- 
cellent asymptotically. This leads us to the conclusion 
that the critical behavior of the 3D EA model is most 
probably universal and that violations of universality as 
they are found in dynamical MC simulation s 12 ' 13 ' 14 ' 16 are 
probably due to difficulties to control the corrections to 
scaling in such simulations. 22 - 23 In Ref. the conclusion 
that the critical behavior of the 3D EA model is univer- 
sal is reached from a high-precision MC study in which 
the static properties of the model with Gaussian and ± J 
couplings are compared. 

We find indications that scaling violations at the given 
bond occupation are in fact reduced in comparison to 
the standard undiluted model, making the study of the 
bond-diluted ± J EA very interesting for various reasons. 

Finally, note that all the figures in this paper where 
the coupling distribution is not explicitly mentioned in 
the caption refer to the data obtained from the study of 
the bond-diluted ± J EA model. 



II. THE MODELS 
A. Bond-diluted ±J EA model 

We mainly consider the 3D bond-diluted EA model 
given by the Hamiltonian 1 



n 



(i) 



where the sum is over all the nearest neighbor pairs of 
a simple cubic 3D lattice of length L. The a x — ±1 are 
Ising spins and the e xy are the bond occupation variables, 



for an empty and e x 



1 for an occupied 



bond. For the couplings J xy we use a binary distribu- 
tion, i.e., J xy — ±1. The fraction of occupied bonds p in 
this study is p — 0.45. We use periodic boundary condi- 
tions to study the Hamiltonian in Eq. ([TJ) . Note that a 
SG transition occurs only at and above a certain critical 
dilution p* that is found to be slightly larger than the 
bond percolation threshold p c — 0.2488 for discrete dis- 
tributions of the couplingsi^i^i This is in contrast to the 
case of continuous coupling distributions, where the SG 
phase appears exactly at p c . A detailed discussion on the 
value of p* will be given in Ref. |23| . 



B. Site-diluted Gaussian EA model 

In order to compare the FSS functions for two differ- 
ent coupling distributions we also consider the 3D site- 
diluted EA model with a Gaussian coupling distribution 
given by the Hamiltonian 1 - 



7~L ^ J Xy£ X £-yt 



(2) 



where the sum is again over all the nearest neighbor pairs 
of a simple cubic 3D lattice of length L. The <j x = ±1 are 
Ising spins and the e x are the site occupation variables, 
i.e., e x — for an empty and e x — 1 for an occupied site. 
For the couplings J xy we use a Gaussian distribution with 
mean zero and variance unity. The fraction of occupied 
sites r in this study is r = 0.625. For this model the SG 
phase is supposed to be present above the site percolation 
threshold r c = 0.3116. We again use periodic boundary 
conditions to study the Hamiltonian in Eq. ^j. 



III. THE FSS METHOD 

Taking a suitably defined finite-volume correlation 
length £(T,L) and a long-range observable 0(T,L) as, 
e.g., X Sg{T 7 L) FSS theory predicts^ 6 . 



Q(T,L) 
0(T, oo) 



= fa [£(T, oo)/L] +OU-",L 



and 



Q(T, sL) 
0(T,L) 



= F o [t(T,L)/L-s]+0{C UJ ,L- 



(3) 



(4) 



where fo & n d Fo are universal FSS functions and s > 1 is 
a scale factor. Eq. ([4]) is an excellent starting point for nu- 
merical investigations of the FSS behavior, as it involves 
only finite-volume quantities taken from a pair of sys- 
tems with sizes L and sL at a given temperature T . The 
knowledge of the universal scaling functions Fo and F% 
allows us to extract critical exponents either by infinite 
volume extrapolatio n 11 ' 20 or by the quotient method i 10 ' 19 
For the quotient method one defines an effective critical 
temperature T* at which the correlation length measured 
in units of the lattice size L are equal for the pair of sys- 
tems, i.e., 



Z{T* c ,L)/L = t{T*,sL)/{ S L). 



(5) 



For an observable O that in the thermodynamic limit 
diverges as t~ x ° (i = T/T c — 1 being the reduced tem- 
perature) the critical exponent xq can be extracted from 
the quotient 



xo/v _ 



0(T*,L) 



+ 0(L- U ) 



(6) 



and hence from Fq. The critical exponent v describing 
the divergence of £ can be, e.g., measured from the FSS 
function of the temperature derivative of £, i.e., from 
Fq t £, which after a short calculation leads to 
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l + —d x F c (x;s)\,+0(L- 



(7) 



with x = £{T,L)/L and x* = t(T*,L)/L. The infinite 
volume extrapolation technique works with data strictly 
above T c and uses the FSS functions and Fo to ob- 
tain the thermodynamic limit of O using an iterative 
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procedure in which the pair £ and O is scaled up from 



sL 



oo as described in detail in 



Ref. l20l . Assuming an ordinary second-order phase tran- 
sition xq can be extracted from a power-law fit to the 
extrapolated data of 0(T, 00). 



IV. ALGORITHM AND COMPUTATIONAL 
DETAILS 

In this MC study we make use of the fact that the 
Hamiltonians in Eqs. ([T]) and @ allow for an efficient 
use of replica cluster algorithms that were originally used 
for the study of the two-dimensional EA model^^ The 
details of the algorithm as they are given below apply 
to the case of the bond-diluted ±J model, whereas a 
discussion of the algorithm in the case of the site-diluted 
Gaussian case can be found in Ref. 0. The basic idea 
of all these algorithms is that two independent replicas 
of the system are simulated simultaneously. Hence each 
lattice site can be associated with two spins, one from 
each replica, and can be in one of the four spin states 

(++), (+-), (-+), and ( ). The (+-) and (-+) sites 

are called A sites and the (++) and ( ) sites B sites. 

On these two different groups of sites it is now possible 
to define clusters in various ways. In the following we 
define the cluster moves used for this simulation. First 
we consider the clusters on the A sites. Following Redner 
et alA we grow the clusters by adding links between two 
A sites with probability p a dd given by 



Padd = 



1 - exp(4:l3E xy ) if E xy < 
else, 



(8) 



where E xy = — J xy e xy a x a y is the energy contribution 
of the given link and (3 is the inverse temperature. We 
repeat the same procedure for the links between the B 
sites. Finally, we flip each of the clusters that are de- 
fined through this procedure with probability 1/2, i.e., 
we make a Swendsen-Wang update. 2 - This choice is influ- 
enced by the findings that even in a diluted ferromagnet 
an update of all clusters can be efficient^ In the pres- 
ence of an external field the algorithm would have to be 
modified slightly as in this case only the clusters on the 
A sites flip freely^ 5 - 

In order to improve the performance of this cluster 
algorithm we also update each replica with parallel tem- 
pering (PT) (Ref. [29|) as proposed by Houdayer^ This 
means that in addition to the two replicas at tempera- 
ture T we introduce replicas on a set of Nt tempera- 
tures (T; with i = 1, . . . , Nt) and allow two replicas at 
Tj and T; + i to exchange their temperatures with the ap- 
propriate probability^ The resulting algorithm is only 
a slight variation of the original replica cluster algo- 
rithm by Swendsen and Wang 4 In contrast to Houdayer's 
algorithm 5 it does not need PT to ensure ergodicity, be- 
cause the different definition of the clusters already pro- 
vides for it. The dilution in the coupling distribution 
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FIG. 1: Qualitative comparison of the integrated autocorrela- 
tion time T au to as a function of temperature between parallel 
tempering (PT) and the replica cluster algorithm (Cluster) 
averaged on the same 20 configurations at L — 10. 



is very important in order to avoid the site percolation 
problem mentioned by Houdayer£ Using, e.g., a nondi- 
luted ±J coupling distribution the cluster move essen- 
tially just swaps the two replicas. That dilution may help 
greatly to speed up cluster algorithms in SG simulations 
has already been noted in the case of the Viana-Bray 
model 30,31 . Finally, we note that it would also be possi- 
ble to define clusters along the lines of Ref. H3, where the 
size of the resulting clusters can be controlled in a more 
general manner. 

An important technical advantage of this algorithm 
with the given discrete distribution of the couplings is 
that it can be completely multispin coded, apart from the 
task where the clusters are identified with the Hoshen- 
Kopelman algorithm. 33 In Fig. [1] we compare the inte- 
grated autocorrelation times of the square of the order 
parameter q 2 (as defined below) of this cluster algorithm 
with the corresponding ones of PT with the same number 
of replicas simulated in total, i.e., with the same number 
of temperature slices, the same set of temperatures, and 
two replicas at each temperature. For this qualitative 
comparison we show an average over the same 20 configu- 
rations on L = 10 lattices. One sees that in the SG phase 
the autocorrelation time of the cluster algorithm for the 
particular choice of temperatures and size is one order 
of magnitude smaller than the one of PT. The compu- 
tational overhead of the cluster move can be kept rather 
small, such that the overall performance remains clearly 
favorable for the cluster algorithm. 

We measure the overlap q x — <j x t x and q = L~ 3 J2 X Qx 
from two replicas a and r with the same couplings 
Jxy^xy We define £(T, L) by the second-moment cor- 
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TABLE I: Simulation parameters for the bond-diluted ±J 
EA model. Number of samples N s , minimal number of 
Monte Carlo sweeps for thermalization Athcrm, the temper- 
ature range T, and the number of temperatures Nt used for 
the PT as a function of the size L. 

where the wave- vector-dependent SG susceptibility is de- 
fined by 

XsG (k) = ^ E E ^[(^x+rMav (10) 
x r 

and where k min = (0, 0, 2n/L) is the smallest non- 
zero wavevector allowed by periodic boundary conditions. 
The SG susceptibility is denned through xsg{T, L) = 
L 3 (q 2 } = xsg(O). We denote the thermal average at 
temperature T by (. . .)t and the average over the disor- 
der realizations by [. . .] av . We also measure the Binder 
ratio defined b y 36 ' 37 

The parameters used for the simulation of the bond- 
diluted ±J EA model are given in Table [H whereas the 
ones used for the simulation of the site-diluted Gaussian 
EA model are given Table |TTJ The number of replicas 
used at each temperature is two. The number of MC 
sweeps for the measurements is equal to the number of 
sweeps used for thermalization. Thermalization has been 
assured by requiring that the observables agree within 
errors in the last three bins of a logarithmic binning. 
The statistical errors of the MC observables are obtained 
from a jackknife estimate i 3 ^ 

V. FSS ANALYSIS AND RESULTS 

In the following we make a detailed FSS analysis of 
the bond-diluted ±J EA model and present results on 
its critical behavior. After that we discuss a few results 
on the site-diluted Gaussian EA model which are mainly 
to motivate that in this model the finite-size effects are 
strong. Its FSS behavior is discussed in more detail in 
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TABLE II: Simulation parameters for the site-diluted Gaus- 
sian EA model. The meaning of the parameters is the same 
as in Table |T] 

Sec. IVI1 where we make a comparison between the FSS 
functions of the two models. An independent determina- 
tion of the critical exponents of the site-diluted Gaussian 
EA model is not in the scope of this paper as the statis- 
tical quality of our data for this model is clearly poorer 
than the one of the data for the bond-diluted ±J EA 
model, as can be seen from Tables HI and ITTl 



A. Bond-diluted ±J EA model 

In Fig. [2] we show the Binder ratio defined in Eq. (fTTj) 
and the correlation length from Eq. (J9j) for the different 
system sizes used in the simulation (see Table [TJ) . The in- 
tersection of the curves for different system sizes defines 
an effective critical temperature that in the infinite vol- 
ume limit converges to T c . The evidence for a SG tran- 
sition is very clear and a relatively precise and robust 
estimate of T c ps 0.66 can be obtained from this data. 
The intersection of the curves is very clean, especially 
also for the Binder ratio, where for other coupling dis- 
tributions the crossing is not very clean i 10 ' 39 Moreover, 
a comparison with the results for the Binder ratio of the 
high-precision study of the undiluted ± J EA model in 
Ref. [lCj gives evidence that the corrections to scaling for 
the bond-diluted ±J distribution used in our study are 
smaller, at least for this observable. The crossing is also 
very clean for the correlation length (apart from the data 
of the smallest lattice L = 4) , but this is less surprising 
as this quantity typically has relatively small corrections 
to scaling*^ 

In Fig.[3]we show the data for the FSS functions and 
F XSG defined in Eq. The FSS Ansatz works very well 
and we find very little differences between the data for 
the different lattice sizes used, indicating that the finite- 
size effects for these functions are already small. The 
data for the pair L = 4 and L = 8 (not shown in Fig. [3]) 
on the other hand show noticeable finite-size corrections. 
In order to extrapolate the data for £ and xsg to infinite 
volume we fit the data of F% and J^sg to the ansatz 
F(x) = 1 + J2i=i n a i ex P( — i/ x ) with n = 5. The fit 
is performed using the data of the largest system sizes, 
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FIG. 2: The correlation length £(T, L) measured in units of 
the lattice size L and the Binder ratio g{T, L) for the different 
system sizes in function of the temperature T. Note that the 
intersection of the different curves is very clean apart from 
that of the smallest system with L — 4 and define T c with 
good precision. The crossings of £(T, L)/L and g(T,L) occur 
very close to each other which is an indication of small scaling 
violations. In order to make the inversion of the ordering 
below and above the phase transition more evident we have 
connected the data points with cubic splines. 



i.e., the data of the pair L = 8 and L = 16. We found 
only slight changes using also the data from the smaller 
systems and the variations we get in the fitted quantities 
are well covered by the errors given in the final results. 
Using the iterative procedure described in Refs. and 
ITT1 we obtain the infinite volume data / XSG which allows 
us to extract the ratio of the critical exponents 7/1/ = 
2 - rj = 2.34 ± 0.02, as asymptotically f a ~ x~ 3io l v 
for x — > 00, if O has a power-law divergence at T c with 
exponent xq (see Fig. [7J. 



Following Ref. [ll| we fit the data extrapolated to infi- 
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FIG. 3: The finite-size scaling functions F x and for the 
scale factor s — 2 from which the critical exponents v and 
7] are obtained. The insets show the same data around the 
effective critical point defined by F% = 2. The cubic splines 
connecting the data points show that the collapse of the data 
for the different sizes is very good. This is a clear sign of small 
scaling violations. 



nite volume to 

£(T) = Cf (T - T D )~ V 

XsG (T) = c x (T-T c )-t 

as well as 

Xsg(0 = ^ 



1 + (T-T C ) 0M 

n 



(12) 



i=0 



l + ^rfWf AW 



i=0 



(13) 



(14) 



with the correction-to-scaling exponents O^ 1 ' and A^^ . In 
particular we have 9^ — wv and A^ ) = lu. We find that 
including a nonleading correction term for the fits of xsg 
(i.e., n = 1 in Eqs. (fill)) and (fT4")0 reduces the dependence 
on the fitting range considerably, whereas for the fits of £ 
the leading correction is enough for obtaining stable and 
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good fits (i.e. 
values 



in Eq. ([12])). We obtain as preferred 



v = 2.17(15), r\ = -0.336(20), 7 = 4.96(30), 
T c = 0.663(6) and lu = 0.7(3) . 

Note that the values we obtain for 7, v and r\ are well 
compatible with the scaling relation 7 = v{2 — 77). The 
value for the correction-to-scaling exponent u> has to be 
taken as an effective value as we do not separate possible 
analytical from nonanalytical corrections to scaling in our 
fits (see Ref. [l8| for a discussion). 

Using on the other hand the quotient method with 
Eqs. ([6]) and ([7]) again for the data of the pair L = 8 and 
L = 16 we obtain 

v = 2.22(15), X] = -0.349(18), 
£(T C ,L)/L = 0.625(10). 

In order to fine tune the ratio £/L on the two lattices 
to be equal and in order to determine the derivative in 
Eq. ([7]) we use a cubic spline fit to the data. The use of 
the cubic spline fits and the remaining scaling violations 
are the main sources of systematic errors. We consider 
them, however, to be smaller than the statistical errors 
in the determination of the universal FSS functions even 
for the rather large number of coupling realizations we 
use in this study. 

Fig. 0] shows the FSS function F g of the Binder ratio 
g. The value of £(T*, L)/L at which F g = 1 gives an al- 
ternative possibility to define an effective critical point. 
The fact that £(T*, L)/L defined from F g = 1 (shown in 
the inset of Fig. 2]) is slightly larger than the one defined 
from F% = 2 (see inset of Fig. [3]) indicates that there are 
remaining corrections to scaling. The value of £(T* ,L)/L 
taken from the pair L = 8 and L = 16 and the condition 
F g = 1 is ^(T*,L)/L = 0.635(15). Comparing to the cor- 
responding value from the pair L = 8 and L = 16 and the 
condition F c = 2, which is £(T*,L)/L = 0.625(10), we 
can conclude that these remaining corrections to scaling 
are small, however. 

In Fig. [5] we show the universal FSS function defined 
by the Binder ratio g(T, L) as a function of the correla- 
tion length £(T,L)/L. This FSS function is particularly 
interesting as these two dimensionless quantities have a 
rather large cross correlation. This means that all the 
data for the sizes from L — 8 to 16 fall on the same curve 
within their statistical errors. The curve provides a one- 
to-one correspondence between the two quantities that 
are often used to determine the critical point in SG sim- 
ulations. Using the value £(T C ,L)/L = 0.625(10) from 
the condition F^ = 2 as our most reliable determination 
of this quantity we obtain a corresponding value of the 
Binder ratio at the critical point of g(T c ) = 0.742(7). 



B. Site-diluted Gaussian EA model 

In Fig. [5] we show the Binder ratio defined in Eq. (fTi"]) 

and the correlation length from Eq. @ for the different 
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FIG. 4: The FSS function of the Binder ratio g(T, L). The 
inset shows the same data around the effective critical point 
defined by F g = 1. The broken lines connecting the data 
points are cubic splines making it again evident that scaling 
violations are small. 
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FIG. 5: The universal FSS function defined by the Binder 
ratio g(T,L) versus the correlation length £(T,L)/L. Note 
that due to the very small scaling violations the plot defines a 
precise mapping between the two quantities that are typically 
used to define the critical temperature. The inset shows the 
same data around the effective critical point. The broken 
line is a polynomial fit through all the data points (i.e., L = 
8, . . . , 16) in the range 0.5 < £(T, L) /L < 0.8. 



system sizes used in the simulation (see Table HI)). In con- 
trast to the data for the bond-diluted ± J model shown 
in Fig. [5] the crossing of the curves for the correlation 
length is no longer very neat. Moreover, the Binder ratio 
does not cross at all for the system sizes and the tem- 
perature range used in this simulation. With increasing 
system sizes the crossing of the correlation length shifts 
from larger T values toward T rs 0.51. Although there 
is no crossing of the Binder ratio we consider the data of 
the correlation length to be clear enough evidence for a 
SG transition around T c w 0.51. We attribute the fact 
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FIG. 6: The same as Fig. [2] but for the site-diluted Gaussian 
coupling distribution. The correlation length shows a cross- 
ing at T c » 0.51. The large shift between the crossing points 
for the smaller system sizes indicates rather large finite-size 
corrections. These large scaling corrections can be seen even 
better for the Binder ratio, as it does not even cross for the 
different system sizes and temperature range used in this sim- 
ulation. The data points are connected by cubic splines. 



that there is no clean intersection to the presence of very 
strong corrections to scaling in this model and we will 
give evidence for this view in Sec. lVIl We think that the 
corrections are so large that nonleading correction terms 
cannot be neglected. 



VI. UNIVERSALITY OF THE FSS FUNCTIONS 

It is obviously interesting to see how the FSS functions 
of these two coupling distributions having such a different 
level of scaling corrections compare to each other and 
that is what we are going to do in this section. But first 
we compare the FSS functions of the bond-diluted ±J 
EA model with the ones of the undiluted ± J EA model. 
In Fig. [7] we compare the data for the FSS functions 
and / XSG defined in Eq. ([3]) obtained in this simula- 
tion with the data from the 3D undiluted ± J EA model 
obtained by Palassini and Caracciolo in Ref. [TJ The 
agreement of the two data sets is excellent giving very 
strong evidence that above the critical dilution p* the 
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FIG. 7: (Color online) A comparison of the universal FSS 
functions f x and f% for the scale factor s = 2 of the bond- 
diluted ±J EA model (circles, L — 4 to L — 16) with the 
corresponding data for the undiluted ±J EA model from 
Palassini and Caracciolo— (triangles, L — 5 to L = 48). The 
insets show the same data in a log- log plot, making evident 
the power-law decay at large £(T, oo)/L. The excellent agree- 
ment between the two data sets gives strong evidence that 
above the critical dilution p* the critical behavior of the 3D 
±J EA model does not depend on bond dilution. 



critical behavior of the ± J EA model does not depend 
on bond dilution. 

In Fig. [8] we compare the data for the FSS function 
defined by the Binder ratio g(T, L) versus the correla- 
tion length £(T,L)/L for the bond-diluted ±J with the 
site-diluted Gaussian coupling distribution. In Fig. [5] we 
showed that this function shows very little scaling correc- 
tions for the bond-diluted ±J EA model. Not surpris- 
ingly this is not the case for the site-diluted Gaussian EA 
model. The curves for increasing system sizes are, how- 
ever, moving closer and closer to the curve defined by the 
data of the bond-diluted ± J EA model. The agreement 
between the curves for the largest few system sizes is very 
good giving support to the claim that the two models fall 
into the same universality class. 

In Fig. [9] we compare the data for the FSS functions 
F(: and F XSG as defined in Eq. (fj| of the bond-diluted 
± J EA model with corresponding ones of the site-diluted 
Gaussian EA model for the pair L = 8 and L = 16. The 
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FIG. 8: (Color online) A comparison of the universal FSS 
function defined by the Binder ratio g(T, L) versus the cor- 
relation length £(T,L)/L (see also Fig. [5jl is given for the 
bond-diluted ±J and the site-diluted Gaussian EA model. 
The figure shows how the data of the site-diluted Gaussian 
model with increasing system size converges toward the curve 
defined by the data of the bond-diluted ± J model. The data 
for the site-diluted Gaussian EA model is connected by cubic 
splines to make the convergence more evident. 



agreement is again very good. Without performing any 
calculation we can read off this figure that the critical 
exponents rj and v of the two models are fully compatible 
with each other within errors. 



VII. DISCUSSION AND CONCLUSIONS 

Replica cluster MC algorithm a 3 ' 4 ' 5 ^ do not only al- 
low for an efficient simulation of the two-dimensional EA 
model, but also for different 3D diluted models with dis- 
order and frustration;- We have applied a variant of such 
an algorithm to the simulation of the 3D bond-diluted EA 
model with binary couplings. The universal FSS func- 
tions are determined with very high statistical accuracy 
and they show little size dependence indicating that the 
corresponding finite-size corrections are becoming small 
when determined from a pair of systems with L = 5 and 
L = 10 or larger. Clearly a change of the behavior of the 
FSS functions at larger sizes than considered here cannot 
completely be excluded, but seems rather improbable. In 
order to obtain good knowledge on the magnitude of the 
remaining scaling corrections a more precise study with 
even more statistics would be needed. Extrapolating the 
data for the correlation length £ and the SG susceptibil- 
ity xsg to infinite volume we obtain the critical expo- 
nents v = 2.17(15), rj = -0.336(20), and the effective 
leading correction-to-scaling exponent uj = 0.7(3). The 
value of 7 = 4.96(30) that we obtain is well compatible 
with values of v, rj and the scaling relation 7 = v(2 — i]) 
connecting these exponents. For the value of the crit- 
ical temperature we get T c — 0.663(6). On the other 



6.0 

Q 5.0 
H 

_$4.0 
g 3 .0 
2.0 
1.0 



2.0 - 



H 

(N 1.5 
H 



1.0 - 



T 1 1 J 1 1 1 J 1 1 1 J 1 1 1 T" 

'I 1 1 1 1 I 1 



5.5 



5.0 



4.5 



A' 



I .... I 







0.60 0.65 





i0 

H — I — I — — I — I — I — — I — I — I — — I — I — I — — h 



• Bond-diluted +/-J 
A Site-diluted Gaussian 



X 











2.05 
2.00 
1.95 
1.90 



(LLl. 



0.55 0.60 0.65 
J 1 1 1 I l_ 



0.0 



0.2 



0.4 0.6 
£(T,L)/L 



0.8 



FIG. 9: A comparison of the finite-size scaling functions F x 
and for the scale factor s — 2 is given for the bond-diluted 
± J and the site-diluted Gaussian EA model for the pair L — 8 
and L = 16. The insets show the same data around the effec- 
tive critical point defined by F% = 2. From the fact that the 
data collapse is very good we can conclude that the estimates 
for the critical exponents v and rj for the two models are the 
same within errors. 



hand using the quotient method we extract v = 2.22(15) 
and 77 = —0.349(18). Other critical quantities we ob- 
tain from this method are £(T C ,L)/L = 0.625(10) and 
using the one-to-one correspondence between the Binder 
ratio and the correlation length (sec Fig. [5]) wc also have 
g(T c ) = 0.742(7). We slightly prefer the values of the 
critical exponents obtained by the quotient method over 
those obtained by the extrapolation to infinite volume 
technique because of two reasons. The first being that 
the fits are clearly more robust using the quotient method 
and the second being that the exponents are extracted 
at the apparent critical point T* which allows for a 
more precise numerical determination of their values i 10 i 19 
Therefore we quote as final values for the critical expo- 
nents of the 3D bond-diluted ±J EA model v = 2.22(15) 
and rj = -0.349(18). 

The comparison of the universal FSS functions and 
fxsc with the ones from the 3D undiluted ± J EA model 
obtained by Palassini and Caracciolo in Ref. [ll| shows 
perfect agreement. This is strong evidence that bond 
dilution is not a relevant perturbation, i.e., it does not 
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change the critical behavior of the standard 3D ± J EA 
model in a large range of dilutions and most probably 
down to the critical dilution p* at which the SG phase 
appears. This conclusion is furthermore supported by the 
excellent agreement between the critical values obtained 
in this study [v = 2.22(15), rj = -0.349(18)] with the 
ones from the high-precision study in Ref. I id, where v = 
2.15(15) and r] = —0.337(15) are given. 

At this point a remark on the difference between the 
critical exponents obtained by Palassini and Caracciolo 11 
is in order, since they give v = 1-8(2) and r\ = —0.26(4) 
using the same infinite volume extrapolation technique 
used in this paper and moreover as they have (within 
errors) the same universal FSS functions and / XSG as 
we have obtained in this paper. We think that this dis- 
crepancy is due to the fact that a precise determination 
of the critical exponents needs data of very high accu- 
racy close to the effective critical point. This was not the 
case in the Palassini and Caracciolo 1 ^ study that involved 
mainly data from the paramagnetic region and therefore 
they were unable to constrain the fits enough to provide 
exact and robust values for the critical exponents. 

The additional comparison of the universal FSS func- 
tions and F XSG of the bond-diluted ± J coupling dis- 
tribution with the corresponding ones of the site-diluted 
Gaussian coupling distribution shows that claims for the 
nonuniversality of the 3D EA model using dynamical 
M C 12 ' 13 i 14 ' 16 are clearly not supported by a careful static 
MC analysis. The agreement of the FSS functions is very 
good and excludes large differences between the critical 
exponents of the two models. Smaller differences at the 
level of 10% obviously cannot be excluded completely 
with the given precision of this study, but looking at the 
complete behavior of the FSS functions (and not only at 
the critical point) the agreement seems too good to leave 
much space for nonuniversal critical behavior of the 3D 
EA model. This is further supported by the agreement 
between the two models regarding the FSS function de- 
fined by the Binder ratio versus the correlation length. 
In the case of the site-diluted Gaussian EA model this 
FSS function clearly is showing large corrections to scal- 
ing for smaller system sizes that, however, asymptotically 
vanish leading most probably to a unique, i.e., universal 
limiting FSS function. We conclude that from static MC 
simulations, especially from the comparison of different 
FSS functions, we get strong evidence for universal crit- 
ical behavior of the SG transition in the 3D EA model. 
The conclusion that the critical behavior of the 3D EA 
model is universal is also reached in Ref. |l7|. There the 
authors perform a high-precision MC simulation to study 



the static properties of the model using mainly a Gaus- 
sian and a ±J distribution for the couplings. It has to 
be noted that the values of the critical exponents they 
find, especially the one for rj of the ±J distribution, is 
not fully consistent with the ones of the present study as 
they find v = 2.44(9) and T) = -0.37(5) for the Gaus- 
sian as well as v = 2.39(5) and r\ = —0.395(17) for the 
± J distribution. We consider this difference to be due to 
the different techniques to extract the critical exponents 
used in their paper, which especially for the value of r/ 
is rather sensitive to nonanalytical as well as analytical 
corrections to scaling. In Ref. |40| the 3D ±J EA model is 
studied on very large lattices with a new technique to ex- 
tract critical exponents. The critical exponents obtained 
in this study are v = 2.72(8) and r\ = -0.4(4) and differ 
again from the ones of our study and also the ones of 
Ref. [lj]. We think that the spread in the above results 
reflects the difficulties to extract critical quantities in the 
3D EA model with high precision and in this respect we 
consider the comparison of the complete FSS functions a 
more reliable and moreover parameter-free tool to check 
for universality. 

Finally, we have noted that the observables we have 
studied in this paper show rather little scaling corrections 
for the 3D bond-diluted ± J EA model at the given bond 
occupation of p — 0.45. This is clearly in accordance with 
the findings of the high-temperature series study in Ref. d 
and the ground state defect energies study in Ref. [9. This 
important physical aspect together with all the technical 
advantages of the replica cluster MC algorithm and the 
possibility of multispin-coding makes this model clearly a 
very interesting candidate for further studies of the static 
properties of spin glasses. 
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